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Suppose a small compact object (black hole or neutron star) of mass m orbits a large black hole 
of mass M 3> m. This system emits gravitational waves (GWs) that have a radiation-reaction effect 
on the particle's motion. EMRIs (extreme-mass-ratio inspirals) of this type will be important GW 
sources for LISA. To fully analyze these GWs, and to detect weaker sources also present in the LISA 
data stream, will require highly accurate EMRI GW templates. 

In this article I outline the "Capra" research program to try to model EMRIs and calculate their 
GWs ab initio, assuming only that m <ti M and that the Einstein equations hold. Because m <^ M 
the timescale for the particle's orbit to shrink is too long for a practical direct numerical integration 
of the Einstein equations, and because this orbit may be deep in the large black hole's strong-field 
region, a post-Newtonian approximation would be inaccurate. Instead, we treat the EMRI spacetime 
' as a perturbation of the large black hole's "background" (Schwarzschild or Kerr) spacetime and use 

^SJ ' the methods of black-hole perturbation theory, expanding in the small parameter m/M. 

The particle's motion can be described either as the result of a radiation-reaction "self-force" 
acting in the background spacetime or as geodesic motion in a perturbed spacetime. Several different 
' lines of reasoning lead to the (same) basic 0{m/M) "MiSaTaQuWa" equations of motion for the 

particle. In particular, the MiSaTaQuWa equations can be derived by modelling the particle as 
either a point particle or a small Schwarzschild black hole. The latter is conceptually elegant, but 
the former is technically much simpler and (surprisingly for a nonlinear field theory such as general 
relativity) still yields correct results. 
^ ' , Modelling the small body as a point particle, its own field is singular along the particle worldline, so 

it's difficult to formulate a meaningful "perturbation" theory or equations of motion there. Detweiler 
and Whiting found an elegant decomposition of the particle's metric perturbation into a singular 
part which is spherically symmetric at the particle and a regular part which is smooth (and non- 
symmetric) at the particle. If we assume that the singular part (being spherically symmetric at the 
particle) exerts no force on the particle, then the MiSaTaQuWa equations follow immediately. 
The MiSaTaQuWa equations involve gradients of a (curved-spacetime) Green function, integrated 
^ ■ over the particle's entire past worldline. These expressions aren't amenable to direct use in practical 

' computations. By carefully analysing the singularity structure of each term in a spherical-harmonic 

expansion of the particle's field, Barack and Ori found that the self-force can be written as an 
00 . infinite sum of modes, each of which can be calculated by (numerically) solving a set of wave 

equations in dimensions, summing the gradients of the resulting fields at the particle position, 
[ and then subtracting certain analytically-calculable "regularization parameters" . This "mode-sum" 

. regularization scheme has been the basis for much further research including explicit numerical 

calculations of the self-force in a variety of situations, initially for Schwarzschild spacetime and 
more recently extending to Kerr spacetime. 

Recently Barack and Golbourn developed an alternative "m-mode" regularization scheme. This 
regularizes the physical metric perturbation by subtracting from it a suitable "puncture function" 
' approximation to the Detweiler- Whiting singular field. The residual is then decomposed into a 

^ , Fourier sum over azimuthal (e""'^) modes, and the resulting equations solved numerically in 2+1 di- 

■ mensions. Vega and Detweiler have developed a related scheme that uses the same puncture-function 

regularization but then solves the regularized perturbation equation numerically in 3-1-1 dimensions, 
avoiding a mode-sum decomposition entirely. A number of research projects are now using these 
puncture-function regularization schemes, particularly for calculations in Kerr spacetime. 

Most Capra research to date has used 1st order perturbation theory, with the particle moving 
on a fixed (usually geodesic) worldline. Much current research is devoted to generalizing this to 
allow the particle worldline to be perturbed by the self-force, and to obtain approximation schemes 
which remain valid over long (EMRI-inspiral) timescales. To obtain the very high accuracies needed 
to fully exploit LISA's observations of the strongest EMRIs, 2nd order perturbation theory will 
probably also be needed; both this and long-time approximations remain frontiers for future Capra 
research. 
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This article is dedicated to the memory of Thomas Radke, my late friend and colleague in many computa- 
tional adventures. 



I. INTRODUCTION 



An EMRI (extreme-mass-ratio inspiral) is a binary black hole (BH) system (or a binary BH/neutron-star system) 
with a highly asymmetric mass ratio. That is, an EMRI consists of a small compact object (a stellar-mass BH or 
neutron star) of mass fiM orbiting a large BH of mass M, with the mass ratio /i <C 1. If the small body were a 
test mass (m = 0), then it would orbit on a geodesic of the large BH. However, if to > 0, then the system emits 
gravitational waves (GWs), and there is a corresponding radiation- reaction influence on the small body's motion. 
Calculating this motion and the emitted GWs is a long-standing research question, and is interesting both as an 
abstract problem in general relativity and as an essential prerequisite for the full success of LISA. LISA is expected 
to observe GWs from many EMRIs with M - IQ^Mq and m ~ lOM© (so that ^J. ~ 10"^) [l|,[2|. To most effectively 
analyze this LISA data - indeed, even to detect much weaker signals that may also be present in the LISA data stream 
- requires accurately modelling the EMRI GWs, particularly the GW phase 

The small body's orbit may be highly relativistic, so post- Newtonian methods (see, for example, 0, section 6.10]; 
and references therein) may not be accurate for this problem. Since the timescale for radiation reaction to shrink 
an EMRI orbit is very long (~ /i~^M) while the required resolution near the small body is very high (~ /iM ), full 
(nonlinear) numerical-relativity methods (see, for example, and references therein) would be both prohibitively 

expensive and insufficiently accurate for this problemjj 

Instead, a variety of other approximation schemes are used to model EMRIs and their GWs. In particular, the 
"Capra" research program^ uses the techniques of BH perturbation theory to model the EMRI spacetime ah initio 
as a perturbation of the massive central BH's Schwarzschild or Kerr spacetime, making no approximations other than 
that the mass ratio <C 1. In particular, the Capra research program doesn't make any slow-motion or weak- field 
approximations. 

In this article I give a relatively non-technical overview of some of the highlights of the Capra research program, 
focusing on those aspects most relevant to explicitly calculating radiation-reaction effects in various physical systems. 
My goal is to give the reader some sense of the "flavor" of Capra research. The reader should have a reasonable 
background in general relativity and, for some parts of sections III Al III Bl and III Dl be familiar with Green-function 
method^H for solving linear partial differential equations (PDEs). The sections of this article are relatively independent 
and, with a few exceptions (which should be obvious from cross-references), can be read in any order. In sections III B I 
and III Dl I have marked certain passages as somewhat more technical (analogous to the "Track 2" of Misner, Thorne, 
and Wheeler [l^); this material may be skipped if the reader so desires. 

This is emphatically not a comprehensive review - there are major areas of the Capra program that I only briefly 
mention, and others which I omit entirely|3 Except for some of the accuracy arguments in section IIVI there's no 
original research in this article. For more detailed and complete information about the Capra research program, the 
reader should consult any of a number of excellent review articles, notably those by Poisson (23l - l25i] Fl Detweiler p6| . 
and Barack [13]. The websites of recent Capra meetings [2^, [1^ also include archives of meeting presentations. 

A key long term goal of the Capra research program is the modelling (and explicit calculation) of highly accurate 
orbital dynamics and GW templates for generic EMRIs. As discussed in section IIVI the highest-accuracy GW 
templates for LISA will probably require carrying BH perturbation theory to at least 2nd order in the mass ratio /z, 
and also using special "long-time" approximation schemes. These are ambitious goals, which are still far from being 
met: most Capra research to date has been devoted to the lesser - but still challenging - problem of trying to model 
strong-fleld EMRI radiation-reaction effects using 1st order perturbation theory and, to the best of my knowledge, no 
Capra GW templates have yet been published. I return to 2nd-order calculations in sections IIVI and |Vl but for the 
rest of this article I consider only Ist-order calculations. 



^ The most asymmetric mass ratio yet simulated with full (nonlinear) numerical relativity is 100:1, i.e., ^ = 10~^ [T^ . A number of 
researchers have attempted to develop special methods to make EMRI numerical-relativity simulations practical, at least for systems 
with "intermediate" mass ratios ^ ~ 10""^. Although promising initial results have been obtained (see, for example, Il7l42ll ). it has not 
(yet) been possible to perform accurate EMRI numerical evolutions lasting for radiation-reaction time scales. 

^ The Capra research program, and the yearly Capra meetings on radiation reaction in general relativity, are named after the late American 
film director Frank Capra, famous for such films as It's a WondcrfuJ Life and Mr. Smith Goes to Washington as well as the World War 
II propaganda series Why We Fight. He owned a ranch near San Diogo and upon his death donated part of this to Caltech. The first 
Capra meeting was held there in 1998. 

^ We say "Bessel function", not "Bessel's function", so logically the reader should be familiar with "Green-function methods", not 
"Green's-function methods. 

I apologise to the reader for any mistakes there may be in this article, and I particularly apologise to anyone whose work I've slighted 
or mischaracterized. I welcome corrections for a future revision of this article. 
^ In particular, Poisson's GR17 plenary lecture [gj contains a short and relatively non-technical review of a large part of the theoretical 
background underlying the Capra research program. I highly recommend this article to the reader seeking somewhat more detail than 
I provide in section [IT] Poisson's lectures [2q| from the 2008 "Mass and Motion" summer school and 11th Capra meeting provides a 
somewhat more detailed presentation of this material, and his Living Reviews in Relativity article |23l gives a lengthy and detailed 
technical account. 
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In almost all Capra calculations to date, the small body is taken to move on a fixed geodesic worldline of the 
background (Schwarzschild or Kerr) spacetime, with radiation- reaction effects being manifest as an 0{fi^) "self- force" 
acting on the small bodyH Alternatively, we can view the small body as moving on a geodesic of a C'(/i)-perturbed 
spacetime. These two perspectives can be shown to be fully equivalent fsoj and are, in some ways, analogous to 
Eulerian versus Lagrangian approaches to fluid dynamics; we can use whichever is more convenient for any given 
calculation. 

Another important choice in self-force analyses is whether to model the small body as a point particle or as a 
nonzero-sized small compact body. Modelling it as a nonzero-sized body is conceptually elegant but technically 
difficult. In contrast, point-particle models are technically simpler but pose difficult conceptual and foundational 
problems. Indeed, in a nonlinear field theory such as general relativity, the very notion of a "point particle" is difficult 
to formulate in a self-consistent manner [3l| . Remarkably, it turns out that these difficulties can be overcome and, 
in fact, point-particle models have been used for the bulk of Capra research to date. I discuss this point further in 
section |lll 

Starting from the Einstein equations, one can derive the generic 0(/i) "MiSaTaQuWa" equations of motion for the 
small body in an arbitrary (strong-field) curved spacetime. These equations give the self-force in terms of a formal 
Green-function integral over the particle's entire past motion and have now been obtained in several different ways, 
using both point-particle and nonzero-sized models of the small body. 

It's usually not possible to explicitly calculate the Green function appearing in the MiSaTaQuWa equations. Instead, 
practical computational schemes are usually based on regularizing the (singular) metric-perturbation equations for 
a point particle; several different ways are now known to do this. The regularized equations can then be solved 
(usually numerically) to actually compute the self-force for a given physical system. Because these calculations are in 
many cases both conceptually difficult and computationally demanding, new techniques are often first developed on 
simpler electromagnetic or scalar-field "model" systems. These retain many of the basic conceptual features of the 
gravitational case while greatly simplifying the gauge choic43 and the resulting computations. 

We can categorize Capra self-force calculations along several dimensions of complexity: 

• The background spacetime may be either Schwarzschild or Kerr. 

• The field equations may be for the scalar-field, electromagnetic, or the full gravitational case. 

• The small body may be stationary, in an equatorial circular orbit, in a generic (non-circular) equatorial orbit, 
or in a fully generic (inclined non-circular) orbit in Kerr spacetime. 

The outline of the remainder of this article is as follows: In section|lI]l discuss some of the key theoretical foundations 
of the Capra program including the Barack-Ori mode-sum regularization, the Detweiler- Whiting decomposition of a 
point particle's metric perturbation, several different derivations of the basic Ist-order "MiSaTaQuWa" equations of 
motion for a small compact body moving in a curved spacetime, the Barack-Golbourn and Vega-Detweiler puncture- 
function regularizations and the self-force computational schemes derived from them, and the decomposition of self- 
force effects into conservative and dissipative parts. In section Hill I summarize a recent self-force calculation of Barack 
and Sago [ssj . which provides an almost complete solution of the Ist-order self- force problem for a particle moving 
on a fixed geodesic orbit in Schwarzschild spacetime. In section ITVl I roughly estimate LISA's accuracy requirements 
for GW templates, and outline some of the issues in trying to model EMRI orbital dynamics for long (orbital-decay) 
times to construct such templates. Finally, in section |V] I summarize the progress of the Capra program to date and 
discuss some of its likely future prospects. 

Throughout this article I use c = G = 1 units and a (— , +, +, +) metric signature. I use the Penrose abstract-index 
notation, with abode as spacetime indices. 5{-) is the Dirac (5-function, r denotes proper time along the small body's 
worldline, and a subscript p denotes evaluation at the small body (particle) 's current position, /i <C 1 is the EMRI 
system's mass ratio and M the central BH's mass. Apart from these, the notation in this article varies somewhat 
from section to section; it's always described at the start of each section. 

II. THEORETICAL BACKGROUND 

In this section I discuss some of the main theoretical background and formalisms which underlie the Capra research 
program^ 



® The small body's mass is O(^), so if it were not constrained to moving on a fixed geodesic worldline, the 0(^^) self-force would give 
rise to an 0{^) "self-acceleration" of the small body away from a geodesic trajectory. 

As discussed by Barack and Ori [33 |. the self-force is highly gauge-dependent in a somewhat unobvious non-tensorial manner. (For 
example, there exists a gauge in which the self-force vanishes. Essentially, the gauge transformation follows the small body as it spirals 
in to the massive BH.) There are thus considerable benefits to computing gauge-invariant effects, an approach particularly championed 
by Detweiler. 

* My exposition in parts of this section draws heavily on that of Poisson's GR17 plenary lecture [24l |. 
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A key early result of Capra research was the derivation in several different ways of the basic Ist-order equations 
of motion for a small compact body moving in a strong-field curved spacetime. These equations were first derived 
in 1997 by Mino, Sasaki, and Tanaka [sJl and Quinn and Wald [sH], and (abbreviating the authors' names) are now 
known as the "MiSaTaQuWa" equations. 

The MiSaTaQuWa equations involve gradients of a curved-spacetime Green function, integrated over the particle's 
entire past worldline. We can rarely calculate the Green function explicitly, so the MiSaTaQuWa equations aren't 
useful for practical computations. In section fll Al I discuss the "mode-sum regularization" computational scheme due 
originally to Barack and Ori [36i] . This scheme regularizes each mode of a spherical- harmonic decomposition of the 
(singular) scalar-field or metric perturbation, then solves numerically for each regularized mode in 1+1 dimensions. 
This scheme has been the basis for much further research, including many practical self-force calculations. 

In 2003 Detweiler and Whiting (37j found a Green-function decomposition - and a corresponding decomposition of 
the metric perturbation due to a small particle - into singular and radiative fields, which greatly aids understanding 
the self- force and related phenomena. In section IIIBI I discuss this decomposition and the related Detweiler- Whiting 
"postulate" concerning the physical significance of the singular and radiative fields. 

In section III CI I outline how the Detweiler- Whiting postulate allows the MiSaTaQuWa equations to be derived via 
modelling the small body as a point particle. This derivation of the MiSaTaQuWa equations is quite simple, but it 
does involve the introduction of point particles and the assumption of the Detweiler- Whiting postulate. 

As an alternative, in section lll DI I outline a different derivation of the MiSaTaQuWa equations, this time modelling 
the small body as a small nonrotating (Schwarzschild) BH. This derivation is technically more difficult than the point- 
particle derivation, but it avoids both the introduction of point particles and the assumption of the Detweiler- Whiting 
postulate. 

Recently researchers have developed several new analyses of the self-force problem, leading to much more satisfactory 
derivations of the MiSaTaQuWa (and analogous) equations. These new analyses are fully rigorous and resolve a number 
of past conceptual difficulties as well as opening promising avenues for further research. I (very) briefly outline these 
analyses in section III El 

Recently two groups have developed alternate "puncture-function" regularization schemes for self-force compu- 
tations. Both schemes first subtract from the physical field a suitable "puncture function" approximation to the 
Detweiler- Whiting singular field, leaving a regular remainder field. Barack, Golbourn, and their coauthors fssU^O] 
developed an "m-mode" regularization scheme which then decomposes the perturbation equation into a Fourier sum 
over azimuthal (e""*^) modes, and finally solves numerically for each mode in 2+1 dimensions. Vega, Detweiler, and 
their coauthors [4ll have developed a different puncture-function regularization scheme that numerically solves 
the regularized equation directly in 3-1-1 dimensions, avoiding a mode-sum decomposition entirely. I describe both of 
these schemes in section III Fl 

The self force can be decomposed into conservative (time-symmetric) and dissipative (time-antisymmetric) parts. 
In section [nG] I discuss this decomposition and its physical significance. 



A. The Barack-Ori Mode-Sum Regularization 

The MiSaTaQuWa equations (discussed further in sections [II Bl III E]) give the self- force in terms of the gradient of a 
curved-spacetime Green function, integrated over the entire past history of the small body. (The integral must be cut 
off inflnitesimally before the small body's current position.) For most physically-interesting systems we can't explicitly 
calculate the Green function, so the MiSaTaQuWa equations aren't useful for practical calculations. Instead, almost 
all practical self-force calculations use other regularized reformulations of the (singular) scalar-field, electromagnetic, 
or metric-perturbation equations. 

Building on earlier suggestions of Ori j4^, ; in 2000 Barack and Ori [3^ proposed a practical "mode-sum" regu- 
larization of the field equations, initially for the model problem of a scalar particle moving in Schwarzschild spacetime. 
Barack and various coauthors [45- 5(],] soon extended this to include electromagnetic and gravitational particles in a 
somewhat wider class of spherically symmetric black hole spacetimes, as well as scalar-field, electromagnetic, and 
gravitational particles in Kerr spacetime [5l|. The mode-sum regularization has been the basis for much further 
research including many practical self-force calculations. 

To explain the mode-sum regularization, I consider the scalar-field case - this contains the essential ideas, but is 
technically much simpler than the full gravitational case. Thus, consider a point particle of scalar charge q, moving 
along a timelike worldline F = F(t) in a background Schwarzschild spacetime with metric gat and covariant derivative 
operator Va. In this section we raise and lower all indices with the background Schwarzschild metric gab, and we take 
(t, r, 6, if) to be the usual Schwarzschild coordinates. 
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We take the scalar field cf) to satisfy the usual scalar wave equation 

□0 = -Anq / — ^ ^ ^ " dr' =: S , (1) 

J^oo V~9 

where the integral extends over the entire worldline of the particle, and we define S to be the source term (right hand 
side). We assume that if the scalar field were regular (non-singular) at the particle, it would exert a force 

Fa^q (Va0)p (2) 

on the particle. 

The scalar wave equation ^ can be formally solved by means of a retarded Green function Q{x,x'), 

/ + 00 
g(x,r(T'))dr' , (3) 
-oo 

where the Green function Q{x^x') satisfies 

Air 

ag{x,x') = -^=d\x^-x"') (4) 

and incorporates the appropriate causality relationships. 

Ignoring some terms which aren't relevant here, the self-force on the particle at the worldline event x = (t, r, 9, ip) 
can then be shown to be given by 

Fa{x)^q^ I VaQ{x-T{T'))dT' , (5) 

J —CO 

where the upper limit x~ means that the integral extends over the entire past worldline of the particle prior to (but 
not including) the event x. 

Now consider a spherical-harmonic decomposition of the scalar field (f) and the self- force Fa^ 

oo +e 

'^(^) = E E ^^™(^' ^)'^'™(*' ^) (6a) 

oo +e 

Fa (^) = E E Y^rn (e, ^)Fi-' {t,r) , (6b) 

and sum over the azimuthal mode number m by defining 

+e 

F^{t, r, 0,ip)= J2 Yern{0, ^)F^"'{t, r) , (7) 

rn — — ^ 

SO that the self-force is given by 

oo 

Fa^Y^^a ■ (8) 

It turns out that each individual spherical-harmonic mode F^ is finite, but the sum over I of these modes in 
diverges. However, by carefully analysing the divergence of the scalar field and its Green function near the particle, 
Barack and Ori showed that the related sum 

oo 

Y,{Fi-[AaL + Ba + CjL]) , (9) 

(where L = ^ + and where the quantities Ba, and Ca are described below) does in fact converge and moreover, 
that the self-force is given by 

oo 

= E (Fa - [AaL + Ba+ Ca/L]) - , (10) 
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where the "regularization parameters" Aa, Ba, Ca, and Za are independent of £ and can be calculated semi-analytically 
as elliptic integrals depending on the worldline F and the worldline event x. 

[More generally, additional even-power terms d'^'^/L'^, D^^L^, d'^'^/L^, . . . , can be added to the AaL + Ba + Ca/L 
inner sum in the mode-sum self-force formula (jlOp (with corresponding adjustments to the definition of Za). As 
discussed by Detweiler, Messaritaki, and Whiting [52| , if the additional regularization parameters Z?!'^' can be explicitly 
calculated (analytically or semi-analytically), then adding these extra terms can greatly accelerate the convergence of 
the sum over £. I discuss the numerical treatment of this infinite sum below.] 

By virtue of the spherical symmetry of the Schwarzschild background, the spherical-harmonic decomposition ^ 
separates the scalar wave equation ([T]). That is, each individual 4>^"^{t,r) now satisfies a linear wave equation in 
1-1-1 dimensions on the Schwarzschild background, 

□0^™ + F,(r)0^" = Se,n{t) S{r - rp{t)) , (11) 

where the potential Vi{r) and source term Sim{t) are known analytically, and where rp(t) is the particle's (known) 
Schwarzschild radial coordinate (which is time-dependent if the worldline is anything other than a circular geodesic 
orbit). The wave equation ()11|) can then be solved numerically in the time domain to find the field 0^™. 

Finally, each individual self- force mode i^f™ can be calculated from the £m component of the basic force law ©, 

- g (Va0'")p , (12) 

and then can be calculated from ([7]). 

Alternatively, we can take a frequency-domain approach, augmenting the spherical-harmonic expansions ([6]) with 
a Fourier transform in time, i.e., we can replace those expansions with 



„ oo +£ 

<^(^) = y E E Y,^{d,^)e'^"^'r'"'{r)dw (13a) 



i=0 m= 

and correspondingly replace the m summation ([7]) with 



^ ^ Y,„,{d,^)e'--'Ff"^{r)diu , (13b) 



r +^ 

F^{t,r,0,^)= / ^ y,„(0,^)e*-'"*F„"^™(r)da. . (13c) 
*^ — 



We also similarly Fourier-transform and spherical-harmonic-expand the source term (right hand side) S in the scalar 
wave equation ([T]), 



oo +£ 



fj2 E ytm{e,ip)e''''"'S^irnir)dio . (13d) 

^ f=0 m=-£ 

The rest of the mode-sum regularization then goes through unchanged for each frequency ui, and each individual 
scalar- field mode (/)'^^™(r) now satisfies the ordinary differential equation 

+ H{r)^^ + K.^(r)r''" = S^i^{r) , (14) 

where the coefficients H{r), V^oiir), and Su,£m{i') are again all known analytically. 

The frequency-domain approach involves much simpler computations (ODEs) than the time-domain approach's 
PDEs. If the particle orbit is circular, then only a single frequency is needed, and the integrals over uj are trivial. On 
the other hand, if the particle orbit is noncircular, many frequencies w may be needed to adequately approximate the 
integrals over uj. 

In any case, once the individual F^™ or F'^^'" are known, there remains the problem that the overall £ sum ([5]) 
is an infinite sum. For computational purposes a finite expression is needed. The solution to this lies in the known 
large-^ behavior [13 of F^^^^ := F^ - [A^L + Ba + Ca/L], 

r,(2) r,(4) r,(6) 

+ ^ + ^ + forlargeL, (15) 
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where the coefficients Di^^ are independent of L. [Of course, if a term D^a'^/ has been added to the AaL-\-Ba + Ca/ L 
inner sum in the mode-sum self- force formula ([TUl) . then that term is absent from the large-L series ([IS]).] We partition 
the overall £ sum ([5]) into a finite "numerical" part and an infinite "tail" part0 

fmax OO 

Fa^Y. <rcg + E ^a.reg , (16) 

where ^max ~ 15-30 is a numerical parameter. Then, once all the ^.^g in the numerical part of the sum are known, 
we can fit some number (typically 2 or 3) of terms in the large-£ series ((T5]) to the numerically-computed i^i^rcg values, 
and use the fitted coefficients {D^''} to estimate the tail term. 

B. The Detweiler- Whiting Decomposition 

Suppose we model the small body as a point particle. (In this section we ignore the point-particle foundational 
issues mentioned in footnote 1141 ) Because the particle's own field is singular along the particle worldline, it's not 
obvious how to write a meaningful "perturbation" theory or formulate equations of motion there. The solution (within 
the framework of modelling the small body as a point particle) is to somehow regularize the field so that we have 
finite quantities to manipulate. 

The regularization of the (retarded) field of a point-particle source near that source is a long-standing problem in 
mathematical physics. Dirac |53j studied this problem for the electromagnetic field of a point charge in flat spacetime 
and found a decomposition of the field into a singular part which is spherically symmetric about the charge, and 
a "radiative" part which is regular at the charge. [We might then assume (postulate) that despite being singular, 
the spherically-symmetric part exerts no force on the charge. I discuss the curved-spacetime generalization of this 



assumption (postulate) below.] Dirac's analysis was generalized to curved spacetime by DeWitt and Brehme [54 1 
(with a correction by Hobbs [5a] 



More recently, Detweiler and Whiting [37[ found a fully satisfactory curved-spacetime decomposition of a point 
particle's field (whether scalar, electromagnetic, or gravitational) into singular and regular parts. Here I describe 
this decomposition for the gravitational case. Suppose we have a point mass m, moving (unaffected by external 
forces apart from self- force effects) with 4- velocity u° along a timelike worldline r = r(r) in a "background" vacuum 
spacetime whose typical radius of curvature in a neighborhood of F is 7?. ^ to. [For a LISA EMRI we might have 
TO ~ IOMq while 7^ ~ IQ^Mq (the massive BH's mass).] 

Suppose ''^^gab is the (vacuum) background metric, i.e., the spacetime metric in the absence of the particle, and Va 
is the corresponding covariant derivative!^ Throughout this section we raise and lower indices with the background 
metric '"'^aft- Let hab be the metric perturbation produced by the particle, so that the physical spacetime metric is 
Qab = ''^hab + hab- We introducc the standard trace-reversed metric perturbation hab = hab — ^^^^gabhc'^ and impose 
the Lorenz gauge condition Vah^^ = on the metric perturbation. 

[The material from this point up to and including the paragraph containing equation (j2ip is somewhat more 
technical than the rest of this article and can be skipped without creating confusion.] 
Up to 0{m/TZ) accuracy, the metric perturbation satisfies the linear (wave) equation 

Wchab+2R\\h,d^Tab , (17) 

where Tab is the particle's ((5-function) stress-energy tensor. 

The (linear) perturbation equation PT|) can be formally solved by introducing a suitable retarded Green function, 
which for events within a normal convex neighborhood of the particle (i.e., events x' which are linked to the particle 
event a; by a unique geodesic) can be written in the Hadamard form 

g'"'c'd'{x,x')=U''\,d>{x,x')6{a{x,x'))+V'''',,d'{x,x')9{-a{x,x')) , (18) 



^ This terminology is perhaps unfortunate: this usage of "tail" has no connection at all to the usage of "tail term" in section III Bl and 
elsewhere in this article. 

The notation (the presence of the prefix on the background metric, but not on the corresponding covariant derivative operator) is 



admittedly somewhat inconsistent here. 
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FIG. 1. This figure shows (in black) the particle's timelike geodesic worldline T, (in red) a nearby "field" event x and its past 

lightcone, (in green) the event x' where this lightcone intersects the worldline F, (in blue) an event x^ on F in the past of x' , 

and (with blue dashed lines) two possible "null then scattering then null" paths from a;^ to a; (these are modelled by the tail 

term (|20|l '). where (heuristically) metric perturbations from the event x'^ scatter off the background curvature before returning 
to event x. 

where 9{-) is the Heaviside step function, Synge's world function a{x,x') is half the squared geodesic distance between 
the events x and x'^^ and where U""^ c'd'{x,x') and V""^ c' d' {x , x') are smooth bitensorso The first term in this 
decomposition is a singular "light-cone part" supported only for ct = 0, i.e., only when x' lies on the past light cone 
of X. The second term is a smooth "tail part" supported only for cr < 0, i.e., only when x' lies within (but not on) 
the past light cone of x. Notice that the tail part has support throughout the entire past history of the particle. 
Heuristically, this is because metric perturbations from any time in the particle's past history can scatter off the 
background spacetime curvature and then return to the event x. 
The resulting formal solution of the perturbation equation ()17p is 



where the S accent (superscript) marks the dummy variable of integration (integrating over the worldline), and where 
the integral is over the particle's entire past history prior to the retarded point x' . Figure [T] illustrates the causal 
relationships between F, x, x' , and x^ . Physically, the tail term ([20]) models the effect of "null then scattering then 
null" paths from x'^ to x (shown in blue in figure [Tl where metric perturbations from the event scatter off the 
background spacetime curvature and then return to the event x. 



As explained very clearly in section 2.1.1 of Poisson's Living Reviews in Relativity article [23l | . Synge's world function a{x,x') has the 
following properties: 

• ct{x, x') < if and only if the geodesic connecting x' to x is timelike, i.e., if and only if x' lies within (but not on) x's past light cone. 

• cr{x, x') = if and only if the geodesic connecting x' to x is null, i.e., if and only if x' lies on x's past light cone. 

• cr{x, x') > if and only if the geodesic connecting x' to x is spacelike, i.e., if and only if x' lies outside x's past light cone. 

For the reader not familiar with bitensors, they are built out of the derivatives da(x, x')/dx and da{x, x')/dx' somewhat analogously to 
the way that "ordinary" tensors are built out of vectors and one-forms. Section 2.1.2 of Poisson's Living Reviews in Relativity article [23| 
gives a brief and very lucid introduction to bitensors. 




■velocity at this retarded 



(19) 




(20) 



singular 



radiative 



FIG. 2. This figure shows the causal properties of the Detweiler- Whiting singular and radiative fields. Each subfigure shows (in 
black) the particle's timelike geodesic worldline F, (in red) a nearby "field" event x and its past and future lightcones, (in blue) 
the events x' (past) and x" (future) where this lightcone intersects the worldline F, and (in green) the region of the worldline 
which affects the Detweiler-Whiting singular or radiative fields. 

The Detweiler-Whiting "singular" part of the metric perturbation is then (defined to be) 

/i°&(x) = — Z^"^,rf,(a;,xMu'='u'^' + — Z^'^^»d»(a;,a;'0u'=\'^''-2TO/ V'', ^(x,T(t''))u'=\'^'' dr'' , (21) 
r r J ^, c^d' 

where (as shown in figure ^ the events x' and x" are respectively the intersections of the past and future light cones 
of the event x with the particle worldline F. 

Detweiler and Whiting showed that the (singular) field hg^ defined in this manner satisfies the same metric- 
perturbation equation ([T7| as the physical (retarded) perturbation ft,"^, and furthermore that h'g' is "just as singular" 
as ft,"* on the worldline F. That is, they showed that the "radiative" field 

/i^''(x) = h''\x) - hf{x) (22) 

is in fact smooth on the worldline F (as well as everywhere else). Notice also that the radiative field satisfies the 
homogeneous form of the metric-perturbation equation ()17|) . 

The Detweiler-Whiting singular and radiative fields have unusual causal properties, illustrated in figure [21 the 
singular field depends on the particle's history only between the events x' and x" , while the radiative field depends 
on the particle's entire past history up to the advanced event x" . However, in the limit that the event x approaches 
the worldline F, the radiative field then depends only on the particle's past history. 

The Detweiler-Whiting singular field h'g'{x) is spherically symmetric at the particle. That is, it can be showrP^ 
that if we average the gradient of this field over a 2-sphere of radius e centered on the particle (as seen in the particle's 
instantaneous rest frame), then take the limit e — > 0, this average vanishes. This motivates the Detweiler-Whiting 
postulate: the singular field exerts no force on the particle; the self-force arises entirely from the action of the (regular) 
radiative field. This postulate gives valuable conceptual insight into how the self force "works". This postulate is 
also closely linked to the MiSaTaQuWa equations (I discuss this in the next paragraph and in section Hi C|) and to 
puncture- function regularizations and computational schemes for the self force (I discuss these in section IIlF[) . 

Unfortunately, because the field is singular at the particle, the simple averaging argument described in the previous 
paragraph doesn't constitute a rigorous proof of the Detweiler-Whiting postulate. However, the Detweiler-Whiting 



See Poisson's Living Reviews in Relativity article | 23] for details. 
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postulate is closely linked to the MiSaTaQuWa equations: if we assume the Detweiler- Whiting postulate, then the 
MiSaTaQuWa equations follow almost immediately via the argument outlined in section III CI Because of this close 
linkage, we can reverse the direction of logical implication and argue that the validity of the MiSaTaQuWa equations 
(which are now well-established via the rigorous derivations discussed in sections fll DI andl ll Ep supports the correctness 
of the Detweiler- Whiting postulate. That is, we can argue that the Detweiler- Whiting postulate must be valid, since 
it is central to a derivation (the one outlined in section III Cp which leads to a correct result (the MiSaTaQuWa 
equations). While not a completely rigorous proof, this argument strongly supports the validity of the Detweiler- 
Whiting postulate. 

Harte |56l - [59j and Pound [60l - l62j have recently given rigorous proofs of the Detweiler- Whiting postulate (somewhat 
generalized in some cases). 



C. Deriving the MiSaTaQuWa Equations via Modelling the Small Body as a Point Particle 

If we ignore the foundational issues of point particles^ then the Detweiler- Whiting postulate provides a relatively 
easy route to the MiSaTaQuWa equations: the (Detweiler- Whiting) statement that the self-force arises solely from 
the action of the (regular) radiative field is equivalent to the statement that the particle moves on a geodesic of the 
background metric perturbed by this radiative field, i.e., '"'f/ah + h^jj- The particle's 4-acceleration is thus given by 

a- ^ ((0)g-^ + u'^u") (iVb/i« - VXc)^'^'' ■ (23) 
It's now easy to show that on the particle's worldline F, 

Vc/l5, - -4TO(w(ai?6)dce + RadbeUc)u'^u' + Vc/l^f . (24) 

Substituting ([M)) into ([^ then gives the MiSaTaQuWa equations 

a- = ((oy^ -t- u'^u") - hlf^)u^u'' . (25) 

where we define /i**^'^ — Vc/i*'^''. The corresponding dynamical equations of motion for the particle are 

u''Vbu" = a" . (26) 

Notice that the metric '•^-'gab + h^b is smooth on the particle worldline F. Moreover, because ^^^gab is a vacuum 
solution of the Einstein equations everywhere and h^j^ satisfies the homogeneous form of the perturbation equation (jl7l) , 
the metric '"-'ffab + h^b ^1^° ^ vacuum solution of the Einstein equations everywhere. This gives what Poisson [2j| 
describes as "a compelling interpretation" to the condition ([23l) : the particle moves on a geodesic of the vacuum 
spacetime with metric '^^^Qab + h^b- Unfortunately, this metric doesn't coincide with the actual physical spacetime 
metric ''^''gab + ^ab- 



D. Deriving the MiSaTaQuWa Equations via Modelling the Small Body as a Black Hole 

This point-particle derivation of the MiSaTaQuWa equations is concise, but depends crucially on the Detweiler- 
Whiting postulate. In this section I outline a different derivation, based on modelling the small body as a BH. 
This derivation doesn't require the assumption of the Detweiler- Whiting postulate but it (this derivation) is techni- 
cally much more involved than the point-particle derivation. This derivation is originally due to Mino, Sasaki, and 
Tanaka [s^l; my presentation here is based on that of Poisson's GR17 plenary lecture |24l | . 

In this section I adopt the same notation as in the point-particle derivation (section III Cp except that the small 
body is no longer modelled as a point particle. Because the small body is (locally) free-falling, its motion is actually 
independent of its internal structure (ignoring spin and tidal effects). This "effacement of internal structure" is a 
fundamental property of general relativity {not shared by most other relativistic gravity theories) and is discussed in 
detail in Damour's fascinating review article in the Three Hundred Years of Gravitation volume In view of this 
property, we are free to choose the small body's internal structure for maximum convenience in our analysis; here we 
choose it to be a nonrotating (Schwarzschild) BH. 



Geroch and Traschen 31] have shown that point particles can not consistently be described by metrics with (5-function curvature tensors. 
More general Colombeau-algebra methods may be able to resolve this problem 63] , but the precise meaning of the phrase "point particle" 
in general relativity remains a very delicate question. 
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Our analysis will be based on matched asymptotic expansions of the spacetime metric: Sufficiently far from the 
small BH (the "far zone"), the metric is that of the background spacetime, perturbed by the presence of the small 
BH, 

5aft = 5'r'^'°""'' + 0(m/r) . (27) 

SufRciently near to the small BH (the "near zone"), the metric is that of the small (Schwarzschild) BH perturbed by 
the tidal field of the background spacetime, 

ffa6=5ff™''^'''''' + 0((r/7e)2) . (28) 

Since m <^ TZ, there exists an intermediate "matching zone" where m/r ^ 1 and r/TZ ^ 1, and hence both the 
expansions (|27p and (|28l) are simultaneously valid. In that region these expansions must represent the same vacuum 
solution of the Einstein equations (modulo gauge choice). The small BH's motion is then determined by the matching 
conditions. 

[The material from this point up to the start of section III D 31 is somewhat more technical than the rest of this 
article and can be skipped without creating confusion.] 

We begin by introducing suitable retarded coordinates centered on the worldline F: w is a backwards null coordinate 
constant on each ingoing null cone centered on F (and is equal to proper time on F), and r is an affine parameter 
on the cone's null generators. In this section I use ijk as Penrose abstract indices ranging over the spatial (non-u) 
coordinates only. The angular coordinates $7* = x'^ /r are constant on each generator (one can think of as spatial 
coordinates on a 2-sphere centered on the small BH). 

1. The Near-Zone Metric 

For present purposes it suffices to approximate the spacetime metric sufficiently near the small BH (the "near zone" ) 
by that of a Schwarzschild BH subject to a quadrupole perturbationl^ In a suitable perturbation {vjX—rCl^) of the 
coordinates (u,x*=rf2*), the null-null component of the near-zone perturbed spacetime metric is 

near ^ _ M _ ^ 



r 



f2 



(^1 - £,,n'n^ + o{{f/nf) , (29) 



where £ij = Cvivj = 0{l/Tl^) are the electric components of the Weyl tensor Cabcd', these components measure the 
tidal distortion induced by the background spacetime. 

2. The Far-Zone Metric 

Sufficiently far from the small BH (i.e., in the "far zone"), the null-null component of the background metric is 

%^^ = -{l + 2ra,^' + r^£,jn'n^) +0{{r/nf) , (30) 

where £ij are once again the electric components of the Weyl tensor, evaluated on F. 

The metric perturbation hat produced by the small BH satisfies the linear perturbation equation (jl7p , with Tab = 
in the far zone. Assuming that the far zone lies within a normal convex neighborhood of the small BH, it can be 
shown that the null-null component of the far-zone perturbed spacetime metric is 

+r{2m£,,n^n^ - 2a,n^ + + Cn"f^O 

+0{{mln){r/nf)+0{{rlnf) . (31) 



This quadrupole term is in general just the leading order in a multipolar expansion. 
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Differentiating the tail term (l20l) . we find tliat in terms of the original (physical) retarded Green function G"'^c'd' {x, x"), 
Kfl is given by 

hit - 4m r V, (g 5 5 - k^^abG' 5,s) ix,x^K^u''^ dr^ , (32) 

J — oo \ a a a / 

where the upper integration limit x' means that the integral extends over the entire past worldline of the small 
BH prior to (but not including) the event x' . By cutting off the integration infinitesimally before x' we include the 
(regular) tail part of the Green function and exclude the (singular) light-cone part. As a result, ft-Jj^',! is finite (although 
generally only C°, i.e., continuous but not differentiable) on the worldline. 

5. Matching 

The coordinate transformation between the near-zone coordinates (w, i*=fri*) and the far-zone coordinates (u, x'^—rCl^) 
can be computed explicitly (up to sufficient orders in the small-in-the-matching-zone quantities m/TZ and r/TZ) in 
terms of h^^^^, its integrals and gradients, and £ij. Using this to transform the far-zone metric into the near-zone 
coordinates gives the null-null metric component 

+o{{m/n){finf)+o{{finf) . (33) 

Requiring that this match the same near-zone metric component (j29p up to 0{m/TZ) now gives the 3- acceleration 
of the small BH's worldline as 

_ l/,tail _ , tail (OA) 

from which the MiSaTaQuWa equations (P5|) follow directly. 

Although the full derivation (including all the steps I've omitted in this brief synopsis) is somewhat lengthy, it can 
be made quite rigorous, requiring no unproven assumptions about the physical system. 

E. Other Derivations of the MiSaTaQuWa Equations 

In this section I briefly mention a number of other derivations of the MiSaTaQuWa equations. In the interests of 
keeping this review both short and broadly accessible, I won't describe any of these derivations in detail. 

As well as the matclied-asymptotic-expansions derivation outlined in section IIIDl Mino, Sasaki, and Tanaka (33 | 
also gave another derivation based on an extension of the electromagnetic radiation-reaction analysis of DeWitt and 
Brehme [si ]. 

Quinn and Wald (35| took an axiomatic approach, showing that the electromagnetic self force can be derived by (i) 
using a "comparison axiom" which relates the electromagnetic force acting on charged particles with the same charge 
and 4-acceleration in two possibly-different spacetimes, and in addition (ii) assuming that in Minkowski spacetime the 
half-advanced, half-retarded electromagnetic field exerts no force on a uniformly accelerating charged particle. Quinn 
and Wald also derived the gravitational self force (the MiSaTaQuWa equations) using a similar set of axioms. 

Building on earlier work by Harte (56l - [59| . Gralla, Harte, and Wald f63| have recently provided a rigorous rederiva- 
tion of the classical (DeWitt-Brehme) electromagnetic self-force based on taking the limit of a 1-parameter family of 
spacetimes corresponding to the small body being "scaled down" in charge and mass simultaneously. (Harte's analysis 
also includes a rigorous proof of a generalized form of the Detweiler- Whiting postulate for the scalar-field and elec- 
tromagnetic cases, and for the linearized Einstein equations.) Gralla and Wald ,|35,] have rederived the gravitational 
self-force (the MiSaTaQuWa equations) based on a similar technique; here the small body is "scaled down" in size 
and mass simultaneously. Both of these derivations are mathematically rigorous and make no assumptions beyond 
the existence and appropriate smoothness and limit properties of the 1-parameter families of spacetimes. 

Pound [60l - [6^ has reviewed various derivations of the MiSaTaQuWa equations (including both the ones I've outlined 
here, and others) and developed several new mathematical techniques for analyzing the self-force problem. Using these. 
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he has rederived the MiSaTaQuWa equation in a highly rigorous manner His analysis includes a rigorous proof 
of the (gravitational) Detweiler- Whiting postulate and also provides many valuable insights into future directions for 
the Capra research program; I outline some of these "future directions" in section IV Bl 

F. Puncture-Function Regularizations 

In this section I describe two recently-developed alternate regularization schemes for self-force computations. Both 
schemes begin by considering a "residual field", defined as the difference between the particle's physical field and 
a suitably chosen "puncture function" which approximates the particle's Detweiler- Whiting singular field near the 
particle. By construction, the residual field is finite (although of limited differentiability) at the particle position, and 
it yields the correct self- force in the force law ([2]). The residual field satisfies a scalar wave equation similar to the 
usual one ((T]), except that by construction the right hand side is now a nonsingular "effective source" that can be 
calculated analytically. I describe the puncture function, the effective source, and the basic outline of how they can 
be used to regularize the (singular) field equation in section IIIF II 

The puncture function and effective source are constructed to have certain specified properties near to the particle. 
Their behavior far from the particle can equivalently be described as either (i) they are undefined far from the particle 
but the computational scheme is formulated so as not to make use of them there, or (ii) they are defined everywhere 
but vanish (or are negligibly small) far from the particle. Following Wardell and his colleagues [66l - l7l| . I use the 
terminology (i); note that some other authors use the terminology (ii). 

Given the puncture function and effective source, Barack and Golbourn's "m-mode" scheme [s^ - lioj does a Fourier 
decomposition of the resulting equation into azimuthal (e*""^) modes, and uses a "world tube" technique to remove 
any dependence on the puncture function or effective source far from the particle. The authors then solve numerically 
for each m-mode of the residual field (using a time-domain numerical evolution in 2-1-1 dimensions for each mode), and 
compute the final self-force by summing over all the modes' contributions. I discuss the Barack-Golbourn m-mode 
scheme further in section IIIF 21 

Vega and Detweiler [41,, ^] take a different approach: Given the puncture function and effective source, they 
introduce a smooth "window function" to remove any dependence on the puncture function or effective source far 
from the particle, then numerically solve the resulting equation directly in 3+1 dimensions. I discuss the Vega- 
Detweiler scheme further in section III F 31 

For either scheme, there are actually many possible choices for the puncture function and effective source. These 
differ in their tradeoffs between the difficulty of analytically calculating the puncture function and effective source, 
and how accurately the puncture function approximates the particle's Detweiler- Whiting singular field (and corre- 
spondingly, how small the effective source is and how smooth the puncture function and effective source are at the 
particle). I discuss this further in section III F 41 

Throughout this section we consider a point particle of scalar charge q, moving along a timelike worldline F in (say) 
Kerr spacetime, whose typical radius of curvature in a neighborhood of F is TZ. 

1. The Basic Puncture-Function Scheme 

In this section I describe the basic puncture-function regularization in its simplest form. This is directly applicable 
to the Barack-Golbourn m-mode scheme discussed in section III F 21 but is slightly modified for the Vega-Detweiler 
scheme discussed in section ITl F 31 

We take the scalar field (j) to satisfy the usual scalar wave equation In general the Detweiler- Whiting singular 
field (f)^ isn't known analytically but, by careful analysis of the scalar field's singularity structure near the particle, 
we can construct approximations to the singular field. Thus, we define an "nth order puncture function" as a 

specific approximation - one that is known analytically - to the Detweiler- Whiting singular field near the particle, 
which satisfies 

^S(«)_^S^0^|^|n-l^ (35) 

near the particle, where A is (roughly) the geodesic distance from the particle (see [40] for a precise definition). Notice 
that at this point, the puncture function need only be defined near the particle; in practice, it's usually only defined 



At the conclusion of his main analysis, Pound writes: "This concludes what might seem to be the most egregiously lengthy derivation 
of the self-force yet performed." . 
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within at most a normal convex neighborhood of the particle worldline. I discuss the construction of the puncture 
function further in section Hi F 41 
We define the "residual" field 

near the particle. The residual field is C"^^ at the particle (and C°° elsewhere near the particle) and satisfies the 
wave equation 

= 5'"'^(") , (37) 

with the "effective source" 5''^'^^"^ given by 

where as in section III Al the integral extends over the entire worldline of the particle. In general the effective source 
is ©(A""^) at the particle (and C°° elsewhere near the particle). 

The subtraction in the effective-source definition ((38)) can't be evaluated numerically (both terms are singular at 
the particle), but it can be evaluated analytically using a (lengthy) series-expansion analysis of the field's singularity 
structure. I discuss this further in section Hi F 41 

If (^-^("^ is a sufficiently good approximation to the true Detweiler- Whiting singular field (j)^ near the particle (i.e., 
if the order n is large enough), and the Detweiler- Whiting postulate holds (i.e., the singular field exerts no force on 
the particle), then it's easy to see that gVc/)^'"' at the particle position gives precisely the desired self-force acting on 
the particle. Thus (apart from the difficulties outlined in the next two paragraphs) the self-force can be calculated by 
analytically calculating the effective source, then numerically solving the wave equations ([37]), and then finally taking 
the gradient of 0^'^"-' at the particle position. 

Accurately solving the wave equation (|57)) is made more difficult by the limited differentiability of the effective 
source and residual field at the particle. With standard finite differencing methods, this limited differentiability limits 
the order of finite-differencing convergence attainable very near the particle. Current research is exploring a variety 
of techniques to alleviate this problem including ignoring it (i.e., simply accepting the lower order of convergence) 
modifying the finite differencing scheme near the particle, and using finite-element or domain-decomposition pseu- 
dospectral methods that naturally accommodate well-localized non-differentiability in the fields [1^ [20|, |42, UM ■ 

Another difficulty with puncture-function regularization schemes is that in general the puncture function and 
effective source are only defined within at most a normal convex neighborhood of the particle whereas the physically 
appropriate boundary conditions for the wave equation (j37|) are applied (to the physical field ((>) at infinity. The 
Barack-Golbourn m-mode scheme and the Vega-Detweiler scheme take very different approaches to resolving this 
difficulty; I discuss these in (respectively) sections fll F 21 and III F 31 below. 



2. The Barack-Golbourn m-mode Scheme 



The Barack-Golbourn m-mode scheme for self-force computation [38|-|4d| begins by defining the puncture func- 
tion and effective source S*^^^^^ exactly as just described (section FlIF 1[) . The authors then decompose the 
residual field and effective source into Fourier series in the azimuthal (ip) direction. 



E ^^^''\t,r,e)e^^^ (39a) 



S''«(^\x)= J2 ^m^"Hi,^,^)e™'^ (39b) 

m— — OO 

Away from the particle, the physical scalar field (f) is smooth and can be similarly decomposed, 

oo 

v)^ (t>mit,r,e)e'^'P . (39c) 



771 — — oo 



^"^ It's not clear to me how much of the overall numerical error occurs within a finite-difference-molecule radius of the particle. If this 
fraction is small at practical grid resolutions, then a lower order of convergence at those few grid points might have only a minor impact 
on the overall numerical accuracy. 
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Since we're working on a Kerr background, the wave equation (|37)) now separates, so that each (complex) residual-field 
m-mode (jhn"'^ (t, r, 0) satisfies a modified wave equation 

□™0f„(") = (40) 

in 2+1 dimensions, where the operator is easily derived analytically and where the effective-source m modes are 
given explicitly by 

S''^^'^\t,r,e) = ^ I S«=^(")(i,r,6',(p')e"'""^' V • (41) 

This integral can be done analytically in some cases, but otherwise must be done numerically. 

There still remains the difficulty that the puncture function and effective source are only defined near the particle, 
while the physical field has has outgoing-wave boundary conditions at infinity. To resolve this problem, the authors 
introduce a worldtube (whose size is a numerical parameter, and shouldn't be "too large" in a sense described below) 
whose interior contains the particle worldline V . The authors then define a new "numerical" field 



,Ar(«) 



6m inside the worldtube 



(42) 



I 0m outside the worldtube , 
and solve numerically for this. The numerical field evidently satisfies the equations 

I outside the worldtube . 
Equivalently, one could say that the authors numerically solve the equations 

□„0«(") = S'°ff(") inside the worldtube (44a) 

Dm 0m = outside the worldtube (44b) 

0^/") = 0m - 0;^(") on the worldtube boundary . (44c) 

The authors solve the piecewise modified wave equation (1401) n umericallv for each m using a standard time-domain 
finite-difference numerical evolution code in 2-t-l dimensions!^ The authors use arbitrary initial data on a large 
domain, in the same manner discussed in section fill Dl below. 

[In a finite-difference numerical code, the piecewise aspect of the equations is trivial to implement [S^H^: the 
code stores the grid function 0m*'"', and for each finite differencing operation, checks if the finite difference molecule 
crosses the worldtube boundary. If so, then the code "adjusts" the grid function values being finite differenced (which 

in this case might well be a temporary copy of a molecule-sized region of the actual grid function 0m*'"'') as appropriate 
using (|44cp .] 

With this scheme neither the puncture function nor the effective source are ever needed more than a short distance 
(the maximum finite-difference molecule size) outside the worldtube. Hence, so long as the worldtube isn't too large, 
it's not a problem that the puncture function and effective source aren't defined far from the particle. Outside the 
worldtube, the piecewise equations (j44l) reduce to □m0m = 0, so it's easy to impose the appropriate outgoing-radiation 
outer boundary conditions. 

Finally, the authors show that the self-force F^^^^ is given by 



oo 



) 

I V ^fYi 



(45) 

r(r) 



where the gradient is evaluated at the particle, and where the real fields 0m'"'* are defined by 



(n) 



(j^m^^ if m = 

2Re (^0m*"'e™^) ff m > . ^^^^ 



Other numerical methods are of course also possible. 
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The infinite sum over m in the self-force law (|45p can be approximated with a finite computation using a tail-fitting 
procedure analogous to that described in section III Al 

The Barack-Golbourn m-mode scheme provides a practical and efficient route to self-force computations for a 
variety of physical systems. It is currently the basis for a number of such calculations. Where the original mode-sum 
scheme described in section III Al reduced the self-force problem to the numerical solution of a 2-dimensional set of 
PDEs in 1-1-1 dimensions^ the m-mode scheme reduces the self-force problem to the solution of a 1-dimensional set 
of PDEs in 2-1-1 dimensions. Both schemes have the major advantage that the problem-domain size, grid resolution, 
and/or other numerical parameters can be varied from one PDE to another. This greatly increases the efficiency of 
the numerical solutions. 



3. The Vega-Detweiler Scheme 

The Vega-Detweiler scheme for self-force computation [4ll, takes a somewhat different approach: it begins by 
defining the puncture function and effective source S''^'^*-"-' exactly as described in section Hi F II The authors 

then introduce a real C°° "window function" W chosen (in a manner described further below) such that 

W ^l + 0{{\/'Rf) (47) 

near the particle, and W ^ Q "sufficiently fast" (i.e., W is either exactly zero or has decayed to a negligible value) 
far from the particle, including in the wave zone and at any event horizon(s) in the spacetime. 

The authors then define the residual field in a manner slightly different from the definition of section Hi F II 
using a subscript W to denote "windowed" quantities, the authors define 

= - VF</)^(") . (48) 
so that the windowed residual field satisfies the wave equation 

^4>T = > (49) 

with the windowed effective source given by 

5f = -□(M^0^(")) - 4^9 r°°^V_£>l (50) 

J-oc V=5 

where once again the integral extends over the entire worldline of the particle. By construction, the residual field and 
effective source so defined have the same continuity properties at the particle as described in section Hi F II 

In the same manner as in section III F II if the puncture function (j)^^'^^ is of sufficiently high order and the Detweiler- 
Whiting postulate holds, then it's easy to see that the windowed residual-field gradient Vc?!)^^"^ at the particle position 
gives precisely the desired self-force acting on the particle. Thus the self-force can be calculated by analytically 
calculating the effective source, then numerically solving the wave equation (1491) in 3-1-1 dimensions with the effective 
source (I50p . then finally taking the gradient of the windowed residual field 0^'"'' at the particle position. 

Since the window function is chosen to approach zero "sufficiently fast" far from the particle, it's not a problem for 
this scheme that the puncture function and effective source aren't defined far from the particle. That is, far from the 
particle we have (either exactly or to an excellent approximation) W = and hence 0^^"'' = (p and S*^^"^ = 0, so the 
wave equation (1491) becomes simply = 0. This also makes it easy to impose the appropriate outgoing-radiation 
outer boundary conditions on cj). 

Like the Barack-Golbourn m-mode scheme, the Vega-Detweiler scheme provides a practical and efficient route to 
self-force computations for a variety of physical systems. The Vega-Detweiler scheme is designed to reduce the self- 
force problem to the numerical solution of a (single) wave equation in 3-1-1 dimensions. This type of problem is quite 
similar to that solved by many existing 3-1-1 numerical relativity codes, so the Vega-Detweiler scheme can often reuse 
existing numerical-relativity codes and/or infrastructure. 



I'm describing the time-domain case here; somewhat similar arguments would also apply to a frequency-domain solution. 
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4- Constructing the Puncture Function and Effective Source 

The key to the success of puncture-function schemes (in either the Barack-Golbourn or Vega-Detweiler variants) 
is the construction of the puncture function This essentially requires a careful local analysis of the field's 

singularity structure near the particle. This can be done exactly only in very simple cases (for example, for a static 
particle in Schwarzschild spacetime). In more general cases, such an analysis uses lengthy series expansions and 
(particularly for higher orders n) is usually done using a symbolic algebra system. Once the puncture function is 
known, the effective source can then be computed (again symbolically) directly from the definition p8|) . The resulting 
algebraic expressions are very lengthy, so usually the symbolic algebra system is also used to directly generate C or 
Fortran for inclusion in a numerical code. 

The difficulty (complexity of the expressions) in computing the puncture function and effective source in this way 
rises very rapidly with the puncture function's order n. In practice, 4th order seems to be both practical and a good 
compromise between the difficulty of computing the puncture function and effective source, the expense of evaluating 
the resulting (machine-generated C or Fortran) expressions in a numerical code, and the differentiability (and hence 
order of accuracy) of the numerical solution. 

Wardell and his colleagues [sgI - ItH have developed efficient software for computing puncture functions and their 
corresponding effective sources at (in theory) any order, and these are now being used in a number of self-force research 
projects. In the interests of brevity I won't try to describe the details of how these puncture functions are calculated, 
but Wardell and Vega [7(| give a very clear description of this. 

G. Conservative versus Dissipative Effects 

The self force can be decomposed into conservative (time-symmetric) and dissipative (time-antisymmetric) parts. 
This decomposition is an important conceptual tool for understanding the physical meaning of the self-force. This 
decomposition is also important for practical computations, for reasons described below. 

To actually compute the conservative and dissipative parts of the self force, consider that thus far, we have used 
solely the retarded scalar field 4>'^°*^ :— (f> (or metric perturbation /i^*^* := hab), and our goal has been to compute the 
corresponding retarded self-force F^^^. :— F^^-^f. If we introduce an advanced scalar field (j>^'^^ (or metric perturbation 
^a'b') ^^'^ the corresponding advanced self-force K^dv (both computed in a manner that's the time-reversal of that for 
the corresponding retarded quantities), then as described in more detail by Dolan and Barack [i^], the conservative 
part of the self-force -F'c'ons and the dissipative part F^^^^ can easily be computed via 

FLs = UPret-F:<i.) ■ (51b) 

This decomposition can also be performed mode-by-mode in a mode-sum or m-mode calculation. In some cases 
there are also ways of computing this decomposition without needing to explicitly compute the advanced self-force; I 
describe one such scheme in section UlI Fl 

The Detweiler- Whiting singular field is time-symmetric, so it cancels out in the subtraction (j51b|) and hence doesn't 
affect the dissipative part of the self-force. This means that the dissipative part can be computed without regularizing 
the singular field, i.e., given a suitable computational scheme, the dissipative part can be computed much more easily 
than the conservative part. In mode-sum and puncture-function regularization schemes, the dissipative part of the 
mode sums also converges much faster (exponentially instead of polynomially) than the conservative part. 

The dissipative part of the self-force directly causes secular drifts in the small bod y's orbital energy, angular mo- 
mentum, and (for non-equatorial orbits in Kerr spacetime) Carter constant. Mino |74l - [78| has argued that dissipative 
self-force alone can be used to calculate the correct long-term (secular) orbital evolution of an EMRI system: the 
conservative part of the self-force appears to cause only quasi-periodic oscillations in the orbital parameters, not long- 
term secular drifts. This "adiabatic approximation" is very useful and can provide a route to EMRI orbital evolution 
that's much simpler and more efficient than the full Capra calculations that are the main subject of this article. 

However, Drasco and Hughes ^] and Pound and Poisson 80, .81] have found that the adiabatic approximation 
isn't as accurate as had previously been thought. In particular, they have found that conservative effects also lead 
to long-term secular changes in the orbital motion. Huerta and Gair [s^l have recently estimated the magnitude 
of these latter effects for a quasicircular EMRI inspiral. In their approximate model of a LISA EMRI whose GWs 
accumulate ~ 10^ radians of phase in the last year of inspiral, conservative effects contribute ^ 20 radians during this 
time interval. Conservative effects are likely to be much larger for eccentric EMRI inspirals[3 While a small fraction 



Building on their self-force calculation for arbitrary bound geodesic orbits in Schwarzschild spacetime |33(| (discussed in detail in 
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of the total phase, even the quasicircular-inspiral conservative effects are still large enough to be easily measurable by 
LISA (in section HV Al I estimate that LISA will be able to detect GW phase differences as small as ^ 10^^ radians). 
Thus conservative effects must be included to model EMRIs sufficiently accurately for LISA. 

Pound and Poisson (s^, HII have also drawn useful distinctions between "adiabatic" , "secular" , and "radiative" 
approximation schemes, which have often been confused in the past. 



III. SELF-FORCE VIA THE BARACK-ORI MODE-SUM REGULARIZATION 

In this section I summarize the recent work of Barack and Sago [s^l in which they calculate the 0{fi) gravitational 
self-force on a particle in an arbitrary (fixed) bound geodesic orbit in Schwarzschild spacetime. 

The calculation is done in the Lorenz gauge, decomposing the metric perturbation due to the particle into tensor 
spherical harmonics, solving for the £ — and £ = 1 harmonics via a frequency-domain method, for the £ > 2 
harmonics by numerically evolving a l-)-l-dimensional wave equation in the time domain, and then computing the 
final self- force using the mode-sum regularization described in section fll Al 

This calculation marks a major milestone in the Capra research program and uses techniques typical of many other 
self-force calculations using time-domain integration of the mode-sum-regularized perturbation equations. This calcu- 
lation also illustrates something of the (high) level of complexity involved in self-force calculations for astrophysically 
"interesting" physical systems - the authors report that (even after many years of preparatory research) it took over 
2 years to develop and debug the techniques and computer code for this calculation. 

In this section I use £m as spherical-harmonic indices, and ijk as "tensor component" indices ranging from 1 to 
10 (these indices are always enclosed in parentheses, and index the individual coordinate components of symmetric 
rank-2 4-tensors). 



A. Particle Orbit 



I take the Schwarzschild line clement to be 

ds^ ^ -f dt^ + f-^ dr^ + r^{d6'^ ^sin^ 6 difi^) (52a) 
^ -f dudv + r'^{de^ +s\T? edi/) , (52b) 

where M is the mass of the Schwarzschild spacetime, / = 1 — 2M/r, (t, r, 6, ip) are the usual Schwarzschild coordinates, 
and (w, u) are null coordinates defined by t; = t + r* and u = t — r,,, where 



2Mlog 



r 

2M ^ 



(53) 



is the "tortise" radial coordinate. 

Without loss of generality the authors take the particle orbit to lie in the equatorial plane 9p — ^. The orbit may 
be parameterized by its (dimensionless) semi-latus rectum p and eccentricity e, defined by 

(54a) 



M(rmin + r, 



max J 

mm 



(54b) 



^min 



are the minimum and maximum r coordinate along the orbit. (For a circular geodesic orbit, 
p = r/M and e = 0.) Figure [3] shows the (p, e) corresponding to unstable, marginally stable, and stable orbits. 

We normalize t (proper time along the particle worldline) to be zero at a (i.e., at some arbitrary) periastron passage 
T = Tyn yn- The particle's geodesic motion x°' = x'^{t;p, e) can then be computed by integrating an appropriate set of 
ODEs [84!|, or semi-analytically in terms of elliptic integrals. 



section [mi l. Barack and Sago |83(| have recently studied conservative self- force effects for eccentric orbits in Schwarzschild spacetime. 
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FIG. 3. This figure shows the (p, e) parameter space for bound geodesic orbits in Schwarzschild spacetime. The region 
p > 6 + 2e where stable geodesic orbits are possible is shown in green. In the region shown in red {p < 6 + 2e) there are no 
stable geodesic orbits, only unstable "plunge" ones. The point at p=6, e=0 marks the innermost stable circular orbit (ISCO). 
The line p = 6 + 2e marks the locus of marginally stable orbits, while zoom- whirl orbits are those just to the right of this line. 



B. Mode-Sum Regularization 

Using the mode-sum regularization discussed in section IlIAl Barack and Ori [s^, Barack and Barack et at. [47| . 
have shown that the (Lorenz-gauge) 4-vector gravitational self-force at any event along the particle's worldline is 
given by 

oo 

^^'^ = E^-g, (55) 
where the "regularized self-force mode" F^^^ is given by 

K'^ = Ffun,± - {^li^ + 5) + ' (56) 

where the "full self-force mode" .Ff'jfj] ^ is computed for each £ as described in section IIII C[ the ± refers to two 
different ways of doing this computation (taking one-sided radial derivatives of the metric perturbation at the particle 
from either the outside or inside), and A^. and B° are "regularization parameters" given semi-analytically in terms 

of certain elliptic integrals of the orbit parameters. Each full self-force mode F^^-^^_^ is itself finite at the particle, but 
their sum diverges whereas the regularized sum (j55p converges. 

The computation of ff°ii ± is based on a tensor-spherical-harmonic decomposition of the Lorenz-gauge metric 
perturbation induced by the particle. Let gab be the background Schwarzschild metric (used to raise and lower 
all indices in this section), g be the determinant of gab, V be the corresponding (background) covariant derivative 
operator, hab be the physical (retarded) metric perturbation due to the particle, h — h^^ be the trace of hab, and 
hab — hab ~ \gabh be the trace-reversed metric perturbation. We take hab to satisfy the Lorenz gauge condition 

V'hab = . (57) 

Let be the particle's 4- velocity. 

To first order in hab, the linearized Einstein equations are then 



y'Vchab + 2R'a'^bhcd = -16nTab , 



(58a) 
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where 



T.. = W (58b) 



'9 

is the particle's ((5-function) stress-energy tensor. 

Due to the symmetry of the Schwarzschild background, the linearized Einstein equations (j58p (together with the 
added constraint-damping terms discussed in section IIII Pp are separable into tensorial spherical harmonics via the 
ansatz 

oo +i 10 

'^^^ = 7E E Y.^^h^'^'-(r,t)Y}l>'-ie,^-r) , (59) 

and similarly for Tat- 

The resulting separated equations take the form of the coupled linear wave equations 



^J^i^)im ^ ^ AAgf + ^ xg^/iO")^™ = Sir - rp) , (60) 

J, a i 

where □ is the 2-D scalar wave operator on the Schwarzschild background, and where A/^^j^, a-nd 5'*^*^^™ are 

given analytically as known functions of the indices (z) and/or (j), position, u", and £ and mW\ 

The authors solve the wave equations ((60)) numerically to obtain the Lorenz-gauge metric perturbation modes 
and their gradients along the particle worldline. I describe this numerical solution in section UlI Dl 



C. The Full Force Modes 



Given the Lorenz-gauge metric perturbation modes in a neighborhood of the particle worldline, the authors 

next compute a set of coefficients /°/)™ defined along the particle worldline in terms of (analytically-known) linear 

combinations of ft.^')^"', 5r±/i^*-'^™, and dt±h^^'^^™' , where the ± corresponds to taking one-sided derivatives from the 
outside or inside of the particle orbit respectively. (Due to the <5- function source term in the wave equation (|60p . 
j^{i)tm typically C° near the particle, i.e., is continuous at the particle worldline but its 1st derivatives have 

a jump discontinuity there.) 

Taking into account the tensor spherical harmonic expansion of hab and Tab (as compared to the scalar spherical 
harmonic expansion implicit in the definition of the f^^y^. ) , the authors then compute 

<l - ^ E f E Kt^'A ^'"(^P- <^.) ' (61) 

P m=-l \p=-3 / 

where each J-^^ is a certain (analytically-known) linear combination of the /"jf^™ with the same £ and m. Because 
of the definition (j6ip - and more generally because of the decomposition of tensor spherical harmonics into scalar 
spherical harmonics - a given full force mode F^^^^^ depends on the Lorenz-gauge metric perturbation modes /i*^*^^ 
for i-i <i' < e+3. 



D. Numerical Solution of the Wave Equations (|6Up 



For each {£, m) the authors solve the 10 coupled wave equations (150]) numerically for the 10 /i^')^™ fields, using 
4th order finite differencing on a uniform characteristic (double-null) (u, u) grid. 

The authors' "diamond integral" finite differencing scheme is adapted from those of [s^, [H]. Since the (S- function 
source term in the wave equation (j60p is nonzero only on the particle's worldline, grid cells away from the worldline 



The reader is warned that my notation here differs from that of Barack and Sago: I make all derivatives explicit in the wave equations II6OI I. 
so that N and M are algebraic coefficients, whereas Barack and Sago use M to denote a single set of Ist-order differential operators 
which contains both the 1st derivative and 0th derivative terms in the wave equations I I6OII . 
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have no source-term contribution, allowing a relatively straightforward finite differencing scheme. The handling of the 
source term for those grid cells which are intersected by the particle worldline - or where the finite difference molecule 
is intersected by the particle worldline - is much more complicated, particularly since (for a non-circular orbit) the 
particle generally crosses grid cells obliquely, with no particular symmetry. 
Several other aspects of the numerical solution are worth of note here: 

• The authors found that a direct numerical solution of the wave equations (|60|) was unstable, with rapidly 
growing violations of the Lorenz gauge constraint ([57)1 . Following Barack and Lousto [13], the authors added 
constraint-damping terms to the evolution equations so as to dynamically damp these gauge violations. 

• The correct initial data for the wave equations (|60|) isn't known. Instead, the authors use zero initial data. This 
results in the evolution initially being dominated by spurious radiation induced by the imperfect initial data. 
Fortunately, this spurious radiation dies out (radiates away) within a few orbital periods, so in a sufficiently 
long evolution its influence eventually becomes negligible^ 

• As always for mode-sum schemes, the numerical calculations are done independently for each i and m. This 
makes the calculation trivial to parallelize. 

• The length of evolution required (or equivalently, given the authors' characteristic grid setup, the size of the 
grid) isn't known a priori. Rather, the evolution must be long enough (the grid must be large enough) for the 
initial-data spurious radiation to have decayed to a sufficiently small level and, more generally, for the /t^*)^"* 
field configuration to have reached an equilibrium. In practice, the authors monitor and its gradient 
along the particle worldline, and stop the evolution once these become periodic (with the particle-orbit period) 
to within a numerical error threshold of ^ lO^"*. If the fields don't meet this criterion before the evolution ends, 
the authors increase the size of the grid and rerun the evolution. 



E. Monopole and Dipole Modes 

The authors were unable to obtain stable numerical evolutions of the wave equations ([SO]) for ^ = or ^ = 1. 
Instead, they (Barack, Ori, and Sago [o^) used a frequency-domain method to solve for /i^*)^™ in these cases. 

Because /i^*)^™ is only C*^ at the particle worldline (/i(*)^™ is continuous at the particle worldline but its 1st derivatives 
have a jump discontinuity there), a naive frequency-domain method would have very poor convergence due to Gibbs- 
phenomenon oscillations. The authors (Barack, Ori, and Sago [s^]) have developed an elegant solution to this problem, 
using the homogeneous modes of the wave equations (j60p as a basis for the numerical solution. They report that 
this "method of extended homogeneous solutions" works very well, with the resulting frequency-domain Fourier sums 
converging exponentially fast to the desired /i^*)^™. 



F. Conservative and Dissipative Parts of the Self- Force 

In general, the technique described in section III Gl for decomposing the self- force into conservative and dissipative 
parts requires explicitly (numerically) computing the advanced metric perturbation h^f'^ as well as the usual retarded 
metric perturbation /ijj"^*. This essentially doubles the overall computational effort. 

As an alternative, the authors describe another way of computing the conservative and dissipative parts of the self- 
force using only the usual (retarded) self-force, assuming only that the particle orbit is periodic with a single intrinsic 
frequency. (This is the case for the authors' system of an arbitrary bound geodesic particle orbit in Schwarzschild 
spacetime, as well as for some types of orbits in Kerr spacetime.) Given this condition, the authors extend an argument 
of Hinderer and Flanagan [9^ to infer that 

(62a) 
(62b) 
(62c) 
(62d) 



^adv(T) = 






+^.'ot(-T) 




+FLi-r) 




-F^A-r) 



Recently Field, Hesthaven, and Lau [88l | suggested that some effects of the spurious radiation would in fact not die out even after 
long evolutions. However, Jaramillo, Sopuerta, and Canizares 89] argue that such "Jost junk solutions" are an artifact of a particular 
(inconsistent) implementation of the (5-function source term. In practice, almost all time-domain mode-sum self-force calculations - 
including the Barack-Sago one being presented here - ignore this issue with no apparent ill effect. Thornburg 90, 91] calculated the 
self-force to < 1 part per million relative accuracy using a time-domain mode-sum code which ignored the possibility of Jost (junk) 
solutions, suggesting that the Jost-solution errors, if present, are very small. 
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Assuming that the usual F^^^ is known for an entire particle orbit, the conservative and dissipative parts of the self- force 
can then be calculated via (ISTj) . As noted in section Hi G[ this decomposition can also be performed mode-by-mode. 



G. The £ Sum 



The mode sum (j55p is an infinite sum. For computational purposes a finite expression is required. To this end, the 
authors partition the mode sum, rewriting (j55p as 

where ^max 15 is a numerical parameter, and compute the first term numerically. (This is the main computation; 
recall that it requires numerically solving the wave equations for < £ < ^max+S.) 

To estimate the second term in the partitioned mode sum the authors consider the conservative and dissipative 
parts of the self- force separately. For the conservative part, the authors make use of the known large-£ asymptotic 
series 

T^a T^a T^a 

pa(l) ^ ^2 _Jd±_ + ^6 ^ . . . .g4^ 

cons.rcg + 1 )2 + + 1 )4 + + 1 )6 + ' 

where the coefficients D)J don't depend on i. The authors least-squares fit the first two terms in this series to the 
numerically computed values of Fc°ons,reg for i^min < ^ < ^maxj where ^min ~ 10 is another numerical parameter. Given 
the fitted coefficients {-Df , -Df }, the second term of the partitioned mode sum (|63|) can then be estimated in terms of 
polygamma functions. 

The dissipative part of the mode sum l|63p converges much faster (in fact, exponentially fast) and is thus much 
easier to handle numerically: in practice, F^^^ieg ^^^^ below the numerical error even before t = ^max, so the second 
term in the partitioned mode sum (|63|) is negligible. 



H. Results and Discussion 



The basic result of the authors' computations is the 4- vector Lorenz-gauge gravitational self-force F° as a function 
of time along (around) the particle orbit. Figure |4] shows an example of these results. 

The authors also studied zoom-whirl orbits, finding and analyzing interesting behavior of the self-force during the 
whirl phase of the orbit. In the interests of brevity, I won't discuss this phenomenon here. 

The authors have also used their code to make the first calculation of the C'(^) self- force corrections to the location 
and angular frequency of the ISCO [o^. This calculation is numerically quite delicate, since the ISCO is a singular 
point in the (p, e) space of particle orbits (figure[3]). The authors used two different techniques to make the calculation, 
and also made a number of other tests to validate the accuracy of their results. They found that self-force effects 
shift the ISCO inwards, from r = 6Af to (6 — 3.269 /i)M, and slightly raise its frequency, from Vl — 1/(6''/^ Af) to 
(1 + 0.4869 n/M)/{6^^'^M). 

These results are of great interest, both as a benchmark of the current state of the Capra research program and 
for comparison with other approaches to modelling EMRI dynamics. Notably, they can usefully be compared with 
post-Newtonian (FN) and effective one-body (EOB) calculations [9^[9^. As well as providing valuable tests of each 
formalism, this can help to "calibrate" various undetermined coefficients in the FN and EOB expansions. 



IV. ACCURACY 



In trying to model EMRI orbital dynamics and calculate EMRI GW templates, it's essential to know what accuracy 
(in the GW phase) is needed, and what accuracy is achievable with various approximation schemes. In this section I 
briefly discuss these issues. 
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FIG. 4. This figure shows the self-force components (dashed line) and _F* (solid line) as a function of scaled Schwarzschild 
time t/M, for a particle orbit with p = 10 and e = 0.5 (so that rmin ~ 6|M and rmax ~ 20M). The orbital period is 434M. 
Notice the slight asymmetry of the self-force with respect to the orbit (for example, F'^ peaks slightly after the particle's 
periastron passage at t/M — 0). This is a genuine physical effect, not a numerical artifact. 



A. The Accuracy Needed by LISA 



Matched filtering of the entire years-long LISA data stream would be impractically expensive for detecting EMRIs 
with hitherto-unknown parameters [93, section 3]. However, once EMRIs have been detected by more economical 
search algorithms (T], Q , precision modelling and matched filtering of the full LISA data sell^ become practical and 
even essential to allow detecting and characterizing weaker sources in the presence of strong EMRIs. 

The strongest LISA EMRIs may have signal/noise ratios of up to p ^ 100 after matched filtering P, so phase 
differences on the order of l/p radians should be just detectable in matched filtering. If we want the maximum possible 
science return from the LISA mission, i.e., if we wish to avoid having this science return limited by the finite accuracy 
of our GW templates, then these templates should have GW phase errors of somewhat less than 10 milliradians. Any 
phase errors larger than this run the risk of significantly increasing the overall parameter-estimation error budget. 
Indeed, if we are lucky and LISA detects a very strong EMRI with (say) p ^ 300, the allowable GW phase errors may 
be even smaller, perhaps ~ 1-2 milliradians. 

As a rough approximation, a typical LISA EMRI accumulates ~ 10^ radians of orbital phase during its last year of 
inspiral [l|, so maintaining a phase error of somewhat less than 10 milliradians implies a fractional error of somewhat 
less than 10 parts per billioro in the instantaneous GW frequency (whose integral gives the cumulative GW phase), 
and hence also in the instantaneous EMRI orbital frequency. 

It's non-trivial to translate "required accuracy in the instantaneous EMRI orbital frequency" into "required accuracy 
in a self-force calculation" , but we can make a crude estimate using the recent analyses of Huerta and Gair [s^, table 1] . 
As discussed in section III Gl for a quasicircular inspiral they estimated that 0{fi^) self- force effects (an ©(/i) ^ 
10~^ fraction of the overall self-force) contribute ~ 20 radians to the cumulative GW phase of our typical LISA 
EMRI. This suggests that a GW phase error tolerance of < 10 milliradians corresponds to a fractional accuracy of 
roughly 

10 milliradians ^ ^ ^ ^ 

20 radians 



Flanagan and Hinderer [9g] have recently found tliat many LISA EMRI inspirals will include several strong transient resonance 
crossings. The EMRI osculating-geodesic orbital state vector exiting such a resonance crossing is a very sensitive function of the 
orbital state vector entering the resonance crossing, with the Jacobian 9(post-resonance state)/9(pre-resonance state) ~ fj.~^^'^ ~ 300 
for a canonical 10^ : IOMq EMRI. This may limit precision modelling and matched filtering to the intervals between strong resonance 
crossings, i.e., to perhaps 0(l/3) of the full EMRI-inspiral data set. 
I use the North American definition that billion = 10^. 
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in the overall self-force. This is only a crude estimate, but it does suggest the general order of magnitude of self-force 
computation accuracy needed to match LISA's data quality for strong EMRIs. 



B. High- Accuracy Capra Computations 

How accurate is (will) (can) a Capra-based GW template be? This obviously depends on many factors. Astrophys- 
ically, there are a variety of possible perturbations to the vacuum-Einstein-e quati ons (Kerr + compact-object) model 
used in most Capra calculations to d ate: magnetic fields, accretion disks [99l4l03j |. or even another supermassive BH 
within a few tenths of a parsec jl04|. It remains an open research problem to incorporate such perturbations within 
Capra models. 

Within vacuum-Einstein-equations Capra models, there are two major sources of error in our models of EMRI 
dynamics and GW emission/propagation: 

• Our Capra computational schemes are based on approximations to the Einstein equations. For example, the use 
of 1st order BH perturbation theory implies fractional errors in self-force effects of at least ©(/x) ^ 10^^ due to 
neglected 2nd-order effects; 2nd order BH perturbation theory should bring these errors downtoO(/x2) ~ 10-^°, 
subject to the long-time-approximation issues discussed in section IIV CI below. No practical 2nd-order Capra 
computational schemes exist yet; I discuss prospects for their construction in the future in section fVBI 

• In practice, we numerically solve our equations using finite-precision arithmetic, and using discrete approxima- 
tions to the ordinary or partial differential equations (ODEs or PDEs). As discussed in section Hi Al frequency- 
domain methods only require integrating ODEs; these are potentially very accurate. For example, Detweiler, 
Messaritaki, and Whiting [s^l achieved fractional numerical errors < 1.5 x 10~^ in a computation of scalar -field 
self-force acting on a scalar particle in a circular orbit in Schwarzschild spacetime, and Blanchet et al. |l05l | 
achieved fractional numerical errors < 10^^'^ in a computation of the gravitational self- force for a point mass in 
a circular orbit in Schwarzschild spacetime. However, time-domain methods require numerically solving PDEs, 
and have typical fractional numerical accuracies of at best < 10^^, although Thornburg [oO, loij achieved frac- 
tional accuracies < 10~^ in a time-domain code by combining adaptive mesh refinement with extended-precision 
floating-point arithmetic. 

For very high accuracy all error sources need to be small, i.e., the EMRI's astrophysical environment must be accu- 
rately modelled and the Capra computational scheme must accurately approximate the Einstein equations (probably 
using 2nd order BH perturbation theory and a long-time approximation of the type described below) and the numerical 
computations must be very accurate. 



C. Long-Time Approximation Schemes 

Most Capra research to date has used 1st order BH perturbation theory and taken the small body to move on a 
fixed (timelike) geodesic worldline in the background Schwarzschild or Kerr spacetime. Assuming some computational 
scheme for the self-force, an obvious improvement is to use the computed self-force to perturb the particle's worldline 
away from being a geodesic, updating a "deviation vector" with the C(/x) equations of motion ([26]). Unfortunately, 
the worldline is fixed by the Ist-order Bianchi identity, so it's not obvious that such a scheme can be self-consistent. 
As very clearly described by Pound ^l| , a "gauge relaxation" technique can be used to allow the particle worldline to 
vary, but such a scheme can still be valid for at most a short time: since the particle's orbit gradually shrinks due to 
GW emission, the particle's orbital phase differs from that of a reference geodesic by ^ 1 radian after the "dephasing 
time", which is quite short - ©(/x"^/^). At this point the deviation vector is large, and the whole approximation 
scheme breaks down. 

A much more sophisticated orbital-evolution scheme is needed to avoid this problem, i.e., to remain accurate for 
the (long) orbital-decay timescales 0{fi~^). Hinderer and Flanagan [ol] and Pound [6l| discuss various aspects of 
how such schemes might be constructed. The detailed definition and implementation of such schemes remains a topic 
for future research. 
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V. SUMMARY AND FUTURE PROSPECTS 
A. Past Light Cone 

A major area of Capra research has long been the effort to analyze the singularity structure of the scalar-field, 
electromagnetic, and gravitational perturbations induced by point charges/masses. Much of our current understanding 
of this structure is based on a decomposition due to Detweiler and Whiting [s^] ■ As discussed in section IIIB| the 
Detweiler- Whiting decomposition splits the perturbation into a singular part (which is, in a suitable sense, spherically 
symmetric at the particle) and a "radiative" part which is finite at the particle. In this context it's common to assume 
the Detwiler- Whiting postulate, which asserts that by virtue of its symmetry the singular field exerts no force on 
the particle; self-force effects arise solely from the particle's interaction with the radiative field. This postulate has 
recently been rigorously proved by Harte (56l - [59| and Pound f60l - [63 [ . 

A variety of lines of reasoning lead to the (same) "MiSaTaQuWa" equations of motion for the Ist-order-perturbation- 
theory self-force acting on a small body moving in an external field. The original derivations of these equations are 
due to Mino, Sasaki, and Tanaka [s^l and Quinn and Wald [sl] (thus the name "MiSaTaQuWa"). More recently, the 
rigorous derivations of Gralla and Wald [65*1, Gralla, Harte, and Wald [gJ], and Pound (60l - [6^ have helped to put the 
MiSaTaQuWa equations on a solid mathematical foundation. 

The notion of "point particle" has serious foundational difficulties in a nonlinear field theory such as general 
relativity jSlJ] (see also footnote [T4l). However, at least in Ist-order perturbation theory these difficulties seem to be 
surmountable. As Poisson writes |23l . section 5.5.4], 

The introduction of a point mass in a nonlinear theory of gravitation would appear at first sight to be 
severely misguided. The lesson learned here is that one can in fact get away with it. The derivation of the 
MiSaTaQuWa equations of motion based on the method of matched asymptotic expansions does indeed 
show that results obtained on the basis of a point-particle description can be reliable, in spite of all their 
questionable aspects. This is a remarkable observation, and one that carries a lot of convenience: It is 
much easier to implement the point-mass description than to perform the matching of two metrics in two 
coordinate systems. 

The MiSaTaQuWa equations involve a curved-spacetime Green function which can only rarely be explicitly cal- 
culated. Instead, almost all practical calculations of self-force effects have returned to the scalar-field. Maxwell, or 
Einstein equations (as appropriate), and regularized the point-particle singularity. 

Barack and Ori [36|] (see also [45l - t50l] ) developed the "mode-sum" regularization, which provides a practical route to 
self-force computations. As discussed in section HL^ this scheme first decomposes the field perturbation into spherical 
harmonics. Each individual spherical-harmonic mode can be calculated by numerically solving a linear wave equation 
in 1+1 dimensions, or by solving an ODE if a frequency-domain approach is used. The self-force is then obtained 
by subtracting certain analytically-calculable regularization parameters from the gradient of each mode's field at the 
particle, and finally summing over all modes. 

The Barack-Ori mode-sum scheme has been the basis for much further research as well as having been used for a 
large number of self-force calculations in various physical systems. In section IIIIl I summarize a noteworthy recent 
self- force calculation using this scheme, due to Barack and Sago [ssj. They calculate the Ist-order-perturbation-theory 
gravitational self-force acting on a particle in an arbitrary bound geodesic orbit in Schwarzschild spacetime, and find 
and analyze a variety of interesting physical effects such as the ©(/i) self- force corrections to the ISCO position and 
orbital frequency. This marks a major milestone in the Capra research program. 

Barack and Golbourn [sl] (see also [H, H^) and Vega and Detweiler [4l| (see also [11]) have recently developed 
"puncture-function" regularization schemes. As discussed in section IIIFl these schemes first subtract a puncture 
function (a suitable analytically-calculable approximation to the Detweiler- Whiting singular field) from the physical 
field near the particle, leaving a finite "residual" field. The Barack-Golbourn "m-mode" scheme decomposes the 
residual field into a Fourier series in the azimuthal {lp) direction, and calculates each azimuthal mode by numerically 
solving a linear wave equation in 2-1-1 dimensions on the Kerr background. The self- force is then obtained by summing 
the field gradient at the particle over all modes. The Vega-Detweiler scheme numerically solves the regularized field 
equation (a linear wave equation) for the residual field in 3-1-1 dimensions, bypassing any mode-sum decomposition. 
Both the Barack-Golbourn m-mode scheme and the Vega-Detweiler scheme are now in use for a variety of self-force 
calculations. 

Because of the complexity of self-force calculations, many techniques have first been developed for scalar-field 
or electromagnetic particles, and then extended to the gravitational case. Historically, self-force calculations first 
considered particles moving in Schwarzschild spacetime, but in recent years a growing number of researchers have 
considered particles in Kerr spacetime. 
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The self force and its effects can usefully be decomposed into conservative (time-symmetric) and dissipative (time- 
asymmetric) parts; the latter are often much easier to calculate. It was once thought that accurate EMRI orbital 
evolutions could be obtained ignoring the conservative part, but Drasco and Hughes [t^ and Pound and Poisson [sollsij l 
have found that this isn't the case. Huerta and Gair 82] have estimated the magnitude of conservative effects as 
^ 20 radians of GW phase (out of a total accumulated GW phase of ^10^ radians) during the final year of a typical 
LISA quasicircular EMRI's inspiral; conservative effects are likely to be much larger for eccentric inspirals. 



B. Future Light Cone 



Many of the topics discussed in this article remain active areas of research. Some of the areas where I expect to see 
major developments are: 

• Further study of the transient resonance crossings recently identified by Flanagan and Hinderer [ol] (described 
in footnotes^. These could have a major impact on LISA EMRI data analysis. 

• Many more self-force calculations for particles orbiting in Kerr spacetime. To this end, several research groups 
are actively working on implementin g and extending the puncture-function computational schemes described in 
section IIIFI Warburton and Barack |l06| have also used the "classic" Barack-Ori mode-sum scheme (discussed 
in section Hi A[) for Kerr calculations. 

• The quantitative comparison of different Capra calculations. Sago, Barack, and Detweiler (soj made a very 
important comparison of this type, showing the consistency of a frequency-domain calculation by Detweiler and 
a time-domain calculation by Sago and Barack, despite these using different gauges. Such comparisons serve as 
valuable checks on all the methods and codes involved. 

• The use of Capra calculations to help det ermine fr ee parameters in post-Newtonian approximation schemes, 
such as the recent work of Blanchet et al. |l05l Il07( . These comparisons also serve as valuable checks on both 
the Capra and post-Newtonian schemes and computations. 

• The development of improved long-time approximation schemes, initially using simple "orbit perturbation" ideas 
as described in section ITVBl and later possibly along the lines suggested by Hinderer and Flanagan [ol] and/or 
Pound [6lj (I briefly discuss the need for these in section HV Cp . 

• The development of analyses, and eventually practical computational schemes, based on 2nd-order BH perturba- 
tion theorv r^PI Pound [6l| has recently reviewed this problem, and has suggested at least one possible route to 
the construction of a practical 2nd-order scheme. The computations required are likely to be very complicated 
(both analytically and numerically) but seem to be possible. 

• The computation of Capra GW templates, and later the incorporation of many of the other developments 
mentioned above into the computation of these templates. 

The result of these and many other developments will (I hope) eventually be the calculation of highly accurate 
EMRI orbital dynamics and GW templates. The numerical calculations involved in doing this will almost certainly 
be very expensive, perhaps comparable or even larger in magnitude to those for a full numerical-relativity simulation 
of a comparable-mass binary BH inspiral/coalescence/ringdown. Thus it won't be practical to calculate in this way 
the huge numbers of GW templates that will be needed for LISA data-analysis template banks. Rather, moderate 
numbers of Capra GW templates will be used to calibrate other (cheaper) approximation schemes (perhaps nth- 
generation descendents of the " kludge" w aveforms used todavFI or perhaps new schemes like the effective one-body 
(EOB) ones described by Yunes |ll5l lll6j ). These cheaper schemes will then be used to generate the actual template 
banks. 

Given the very talented people working on these problems, I predict that Capra EMRI GW templates meeting the 
accuracy goal described in section IIVBI (< 10 milliradians of phase error over a full million-radian inspiral) will be 
published within 10 years of this article's appearance. I hope many readers of this article will participate in this effort 
and that within most of our lifetimes we will see actual LISA data being filtered with these templates. There is much 
work to do. 



Rosenthal [lOSHlllll has obtained a formal expression for the 2nd-order self-force, but unfortunately this is in a gauge which is very 
inconvenient for practical calculations (in this gauge the Ist-order self- force vanishes). 

As well as their use for EMRIs, 2nd-order schemes would be of great value in modelling intermedia te mass r atio inspirals. 
For an introduction to kludge waveforms see, for example, Barack and Cutler or Babak et al. |ll3l . lll3| . 
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